function least_squares_applications()
% LEAST_SQUARES_APPLICATIONS 最小二乘应用实例
% 
% 展示QR分解在各种最小二乘问题中的应用

fprintf('=== 最小二乘应用实例 ===\n\n');

%% 1. 多项式拟合
fprintf('1. 多项式拟合\n');
fprintf('目标：使用QR分解进行多项式拟合\n\n');

% 生成带噪声的多项式数据
rng(42);
x_data = linspace(-2, 2, 50)';
true_coeffs = [1, -0.5, 2, 0.3];  % 真实系数 [常数项, x, x^2, x^3]
y_true = polyval(fliplr(true_coeffs), x_data);
noise_level = 0.2;
y_data = y_true + noise_level * randn(size(y_true));

fprintf('数据点数: %d\n', length(x_data));
fprintf('真实多项式: f(x) = %.1f + %.1fx + %.1fx² + %.1fx³\n', true_coeffs);
fprintf('噪声水平: %.1f\n', noise_level);

% 构造Vandermonde矩阵
degrees = [1, 2, 3, 4, 5];
fprintf('\n不同阶数多项式拟合结果:\n');
fprintf('阶数    条件数        残差        拟合系数\n');
fprintf('----    ------        ----        --------\n');

for degree = degrees
    % 构造设计矩阵
    A_poly = zeros(length(x_data), degree + 1);
    for i = 0:degree
        A_poly(:, i+1) = x_data.^i;
    end
    
    % QR分解求解
    [Q_poly, R_poly] = qr(A_poly, 0);
    coeffs_qr = R_poly \ (Q_poly' * y_data);
    
    % 计算残差
    y_fit = A_poly * coeffs_qr;
    residual = norm(y_data - y_fit);
    cond_num = cond(A_poly);
    
    fprintf('%4d    %8.2e    %8.4f    [', degree, cond_num, residual);
    fprintf('%.3f ', coeffs_qr);
    fprintf(']\n');
    
    % 对于3阶多项式，比较不同方法
    if degree == 3
        fprintf('\n3阶多项式不同求解方法比较:\n');
        
        % 正规方程
        coeffs_normal = (A_poly' * A_poly) \ (A_poly' * y_data);
        residual_normal = norm(y_data - A_poly * coeffs_normal);
        
        % SVD方法
        [U_svd, S_svd, V_svd] = svd(A_poly, 0);
        coeffs_svd = V_svd * (S_svd \ (U_svd' * y_data));
        residual_svd = norm(y_data - A_poly * coeffs_svd);
        
        fprintf('方法        残差        系数差异\n');
        fprintf('--------    --------    --------\n');
        fprintf('QR分解      %8.4f    基准\n', residual);
        fprintf('正规方程    %8.4f    %8.2e\n', residual_normal, norm(coeffs_qr - coeffs_normal));
        fprintf('SVD方法     %8.4f    %8.2e\n', residual_svd, norm(coeffs_qr - coeffs_svd));
    end
end

%% 2. 线性回归
fprintf('\n2. 多元线性回归\n');
fprintf('目标：使用QR分解进行多元线性回归分析\n\n');

% 生成回归数据
n_samples = 100;
n_features = 5;
rng(123);

% 真实回归系数
beta_true = [2, -1.5, 0.8, -0.3, 1.2]';
intercept_true = 0.5;

% 生成特征矩阵
X_reg = randn(n_samples, n_features);
% 添加一些相关性
X_reg(:, 2) = 0.7 * X_reg(:, 1) + 0.3 * randn(n_samples, 1);

% 生成响应变量
y_reg = intercept_true + X_reg * beta_true + 0.3 * randn(n_samples, 1);

% 添加截距项
X_reg_with_intercept = [ones(n_samples, 1), X_reg];
true_params = [intercept_true; beta_true];

fprintf('回归问题设置:\n');
fprintf('样本数: %d\n', n_samples);
fprintf('特征数: %d (加截距项)\n', n_features + 1);
fprintf('设计矩阵条件数: %.2e\n', cond(X_reg_with_intercept));

% QR分解求解
[Q_reg, R_reg] = qr(X_reg_with_intercept, 0);
beta_hat = R_reg \ (Q_reg' * y_reg);

% 计算统计量
y_pred = X_reg_with_intercept * beta_hat;
residuals = y_reg - y_pred;
rss = sum(residuals.^2);  % 残差平方和
tss = sum((y_reg - mean(y_reg)).^2);  % 总平方和
r_squared = 1 - rss / tss;  % 决定系数

fprintf('\n回归结果:\n');
fprintf('参数        真实值    估计值    误差\n');
fprintf('--------    ------    ------    ----\n');
param_names = {'截距', '特征1', '特征2', '特征3', '特征4', '特征5'};
for i = 1:length(beta_hat)
    error_i = abs(true_params(i) - beta_hat(i));
    fprintf('%-8s    %6.2f    %6.2f    %4.2f\n', param_names{i}, true_params(i), beta_hat(i), error_i);
end

fprintf('\n模型评估:\n');
fprintf('R²: %.4f\n', r_squared);
fprintf('残差标准误: %.4f\n', sqrt(rss / (n_samples - n_features - 1)));
fprintf('F统计量: %.2f\n', (r_squared / n_features) / ((1 - r_squared) / (n_samples - n_features - 1)));

%% 3. 非线性最小二乘（线性化）
fprintf('\n3. 非线性最小二乘（线性化方法）\n');
fprintf('目标：拟合指数衰减模型\n\n');

% 生成指数衰减数据: y = a * exp(-b * x) + c
t_exp = linspace(0, 5, 30)';
a_true = 3.0;
b_true = 0.8;
c_true = 0.5;
y_exp_true = a_true * exp(-b_true * t_exp) + c_true;
y_exp_data = y_exp_true + 0.1 * randn(size(y_exp_true));

fprintf('指数衰减模型: y = a*exp(-b*t) + c\n');
fprintf('真实参数: a=%.1f, b=%.1f, c=%.1f\n', a_true, b_true, c_true);

% 线性化方法1：对数变换（假设c=0）
y_exp_shifted = y_exp_data - min(y_exp_data) + 0.1;  % 确保正值
log_y = log(y_exp_shifted);
A_log = [ones(length(t_exp), 1), -t_exp];

[Q_log, R_log] = qr(A_log, 0);
params_log = R_log \ (Q_log' * log_y);
a_est_log = exp(params_log(1));
b_est_log = params_log(2);

fprintf('\n对数线性化结果:\n');
fprintf('估计参数: a=%.2f, b=%.2f\n', a_est_log, b_est_log);

% 线性化方法2：差分方法
% 使用Gauss-Newton方法的一步近似
% 初始猜测
a_init = max(y_exp_data);
b_init = 0.5;
c_init = min(y_exp_data);

% 构造雅可比矩阵
J_exp = zeros(length(t_exp), 3);
for i = 1:length(t_exp)
    exp_term = exp(-b_init * t_exp(i));
    J_exp(i, 1) = exp_term;                    % ∂f/∂a
    J_exp(i, 2) = -a_init * t_exp(i) * exp_term;  % ∂f/∂b
    J_exp(i, 3) = 1;                           % ∂f/∂c
end

% 残差
r_exp = y_exp_data - (a_init * exp(-b_init * t_exp) + c_init);

% Gauss-Newton步
[Q_gn, R_gn] = qr(J_exp, 0);
delta_params = R_gn \ (Q_gn' * r_exp);

a_est_gn = a_init + delta_params(1);
b_est_gn = b_init + delta_params(2);
c_est_gn = c_init + delta_params(3);

fprintf('\nGauss-Newton一步结果:\n');
fprintf('估计参数: a=%.2f, b=%.2f, c=%.2f\n', a_est_gn, b_est_gn, c_est_gn);

%% 4. 约束最小二乘
fprintf('\n4. 约束最小二乘\n');
fprintf('目标：求解带等式约束的最小二乘问题\n\n');

% 问题: min ||Ax - b||² subject to Cx = d
m_const = 20;
n_const = 8;
p_const = 3;  % 约束数量

rng(456);
A_const = randn(m_const, n_const);
x_true_const = randn(n_const, 1);
b_const = A_const * x_true_const + 0.1 * randn(m_const, 1);

% 等式约束
C_const = randn(p_const, n_const);
d_const = C_const * x_true_const;  % 确保真实解满足约束

fprintf('约束最小二乘问题:\n');
fprintf('变量数: %d\n', n_const);
fprintf('观测数: %d\n', m_const);
fprintf('约束数: %d\n', p_const);
fprintf('约束矩阵条件数: %.2e\n', cond(C_const));

% 方法1：拉格朗日乘数法
% 构造KKT系统: [A'A  C'; C  0] [x; λ] = [A'b; d]
KKT_matrix = [A_const' * A_const, C_const'; C_const, zeros(p_const, p_const)];
KKT_rhs = [A_const' * b_const; d_const];

[Q_kkt, R_kkt] = qr(KKT_matrix);
solution_kkt = R_kkt \ (Q_kkt' * KKT_rhs);
x_constrained = solution_kkt(1:n_const);
lambda_mult = solution_kkt(n_const+1:end);

% 验证约束满足
constraint_violation = norm(C_const * x_constrained - d_const);
objective_value = norm(A_const * x_constrained - b_const)^2;

fprintf('\n拉格朗日乘数法结果:\n');
fprintf('约束违反: %.2e\n', constraint_violation);
fprintf('目标函数值: %.4f\n', objective_value);
fprintf('拉格朗日乘数: [');
fprintf('%.3f ', lambda_mult);
fprintf(']\n');

% 方法2：零空间方法
% 找到C的零空间
[Q_null, R_null] = qr(C_const');
Z_null = Q_null(:, p_const+1:end);  % 零空间基

% 特解：满足约束的解
x_particular = pinv(C_const) * d_const;

% 在零空间中求解
A_null = A_const * Z_null;
b_null = b_const - A_const * x_particular;

[Q_null_ls, R_null_ls] = qr(A_null, 0);
y_null = R_null_ls \ (Q_null_ls' * b_null);

x_null_space = x_particular + Z_null * y_null;

fprintf('\n零空间方法结果:\n');
fprintf('约束违反: %.2e\n', norm(C_const * x_null_space - d_const));
fprintf('目标函数值: %.4f\n', norm(A_const * x_null_space - b_const)^2);
fprintf('解的差异: %.2e\n', norm(x_constrained - x_null_space));

%% 5. 加权最小二乘
fprintf('\n5. 加权最小二乘\n');
fprintf('目标：处理异方差数据\n\n');

% 生成异方差数据
n_weighted = 50;
x_weight = linspace(0, 10, n_weighted)';
y_weight_true = 2 + 1.5 * x_weight - 0.1 * x_weight.^2;

% 异方差噪声：方差随x增加
sigma_weight = 0.1 + 0.05 * x_weight;
y_weight_data = y_weight_true + sigma_weight .* randn(n_weighted, 1);

% 权重矩阵（方差的倒数）
W = diag(1 ./ sigma_weight.^2);

fprintf('加权最小二乘设置:\n');
fprintf('数据点数: %d\n', n_weighted);
fprintf('权重范围: [%.2f, %.2f]\n', min(diag(W)), max(diag(W)));

% 设计矩阵（二次多项式）
A_weight = [ones(n_weighted, 1), x_weight, x_weight.^2];

% 普通最小二乘
[Q_ols, R_ols] = qr(A_weight, 0);
beta_ols = R_ols \ (Q_ols' * y_weight_data);
residual_ols = norm(y_weight_data - A_weight * beta_ols);

% 加权最小二乘
% 变换: W^(1/2) * A * beta = W^(1/2) * y
W_sqrt = sqrt(W);
A_weighted = W_sqrt * A_weight;
y_weighted = W_sqrt * y_weight_data;

[Q_wls, R_wls] = qr(A_weighted, 0);
beta_wls = R_wls \ (Q_wls' * y_weighted);
residual_wls = norm(y_weighted - A_weighted * beta_wls);

fprintf('\n拟合结果比较:\n');
fprintf('方法        残差        系数 [常数, 线性, 二次]\n');
fprintf('--------    --------    -------------------------\n');
fprintf('普通LS      %8.4f    [%6.3f, %6.3f, %6.3f]\n', residual_ols, beta_ols);
fprintf('加权LS      %8.4f    [%6.3f, %6.3f, %6.3f]\n', residual_wls, beta_wls);
fprintf('真实值      %8s    [%6.3f, %6.3f, %6.3f]\n', '-', 2.0, 1.5, -0.1);

%% 6. 总体最小二乘
fprintf('\n6. 总体最小二乘\n');
fprintf('目标：处理设计矩阵也有误差的情况\n\n');

% 生成数据（设计矩阵和响应都有误差）
n_tls = 30;
x_tls_true = linspace(1, 5, n_tls)';
slope_true = 2.5;
intercept_true = 1.0;
y_tls_true = intercept_true + slope_true * x_tls_true;

% 添加误差
x_tls_obs = x_tls_true + 0.1 * randn(n_tls, 1);
y_tls_obs = y_tls_true + 0.15 * randn(n_tls, 1);

fprintf('总体最小二乘设置:\n');
fprintf('数据点数: %d\n', n_tls);
fprintf('真实参数: 截距=%.1f, 斜率=%.1f\n', intercept_true, slope_true);
fprintf('x误差标准差: 0.1\n');
fprintf('y误差标准差: 0.15\n');

% 普通最小二乘（忽略x的误差）
A_tls_ols = [ones(n_tls, 1), x_tls_obs];
[Q_tls_ols, R_tls_ols] = qr(A_tls_ols, 0);
params_ols_tls = R_tls_ols \ (Q_tls_ols' * y_tls_obs);

% 总体最小二乘（SVD方法）
% 构造增广矩阵 [X, -y]
augmented_matrix = [x_tls_obs, -y_tls_obs];
[U_tls, S_tls, V_tls] = svd(augmented_matrix);

% 最小奇异值对应的右奇异向量给出解
v_min = V_tls(:, end);
slope_tls = v_min(1) / v_min(2);

% 计算截距
intercept_tls = mean(y_tls_obs) - slope_tls * mean(x_tls_obs);

fprintf('\n拟合结果比较:\n');
fprintf('方法        截距      斜率      截距误差    斜率误差\n');
fprintf('--------    ------    ------    --------    --------\n');
fprintf('普通LS      %6.3f    %6.3f    %8.3f    %8.3f\n', ...
        params_ols_tls(1), params_ols_tls(2), ...
        abs(params_ols_tls(1) - intercept_true), abs(params_ols_tls(2) - slope_true));
fprintf('总体LS      %6.3f    %6.3f    %8.3f    %8.3f\n', ...
        intercept_tls, slope_tls, ...
        abs(intercept_tls - intercept_true), abs(slope_tls - slope_true));

%% 7. 递归最小二乘
fprintf('\n7. 递归最小二乘\n');
fprintf('目标：在线更新最小二乘估计\n\n');

% 生成流数据
n_total = 100;
n_features_rls = 3;
rng(789);

% 真实参数（可能时变）
beta_true_rls = [1.5, -0.8, 2.2]';

% 生成批量数据用于比较
X_batch = randn(n_total, n_features_rls);
y_batch = X_batch * beta_true_rls + 0.2 * randn(n_total, 1);

fprintf('递归最小二乘设置:\n');
fprintf('总数据点数: %d\n', n_total);
fprintf('特征数: %d\n', n_features_rls);
fprintf('真实参数: [%.1f, %.1f, %.1f]\n', beta_true_rls);

% 批量最小二乘（基准）
[Q_batch_rls, R_batch_rls] = qr(X_batch, 0);
beta_batch = R_batch_rls \ (Q_batch_rls' * y_batch);

% 递归最小二乘
beta_rls = zeros(n_features_rls, 1);  % 初始估计
P_rls = 1000 * eye(n_features_rls);   % 初始协方差矩阵

fprintf('\n递归更新过程 (每10步显示):\n');
fprintf('步数    当前估计 [β1, β2, β3]        估计误差\n');
fprintf('----    -------------------------    --------\n');

for k = 1:n_total
    x_k = X_batch(k, :)';
    y_k = y_batch(k);
    
    % RLS更新公式
    gain = P_rls * x_k / (1 + x_k' * P_rls * x_k);
    beta_rls = beta_rls + gain * (y_k - x_k' * beta_rls);
    P_rls = P_rls - gain * x_k' * P_rls;
    
    if mod(k, 10) == 0 || k <= 5
        estimation_error = norm(beta_rls - beta_true_rls);
        fprintf('%4d    [%6.3f, %6.3f, %6.3f]    %8.4f\n', ...
                k, beta_rls, estimation_error);
    end
end

fprintf('\n最终结果比较:\n');
fprintf('方法        参数估计 [β1, β2, β3]        估计误差\n');
fprintf('--------    -------------------------    --------\n');
fprintf('批量LS      [%6.3f, %6.3f, %6.3f]    %8.4f\n', ...
        beta_batch, norm(beta_batch - beta_true_rls));
fprintf('递归LS      [%6.3f, %6.3f, %6.3f]    %8.4f\n', ...
        beta_rls, norm(beta_rls - beta_true_rls));

%% 总结
fprintf('\n=== 最小二乘应用总结 ===\n');
fprintf('1. 多项式拟合: QR分解提供数值稳定的解\n');
fprintf('2. 线性回归: 标准的统计建模工具\n');
fprintf('3. 非线性最小二乘: 线性化方法的基础\n');
fprintf('4. 约束最小二乘: 处理等式约束的优化问题\n');
fprintf('5. 加权最小二乘: 处理异方差数据\n');
fprintf('6. 总体最小二乘: 考虑设计矩阵误差\n');
fprintf('7. 递归最小二乘: 在线学习和自适应估计\n\n');

fprintf('QR分解在最小二乘中的优势:\n');
fprintf('• 数值稳定性优于正规方程\n');
fprintf('• 自然处理秩亏问题\n');
fprintf('• 便于增量更新\n');
fprintf('• 提供几何直观\n');
fprintf('• 适合各种约束和扩展\n\n');

end